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[i]  Quasi-linear  diffusion  by  cyclotron-resonant  plasma  waves  is  likely  a  key  ingredient 
of  the  behavior  of  electrons  in  the  Earth’s  radiation  belts.  Multidimensional  dynamical 
simulations  are  under  development,  which  require  the  diffusion  coefficients  to  be 
evaluated  quickly  as  well  as  accurately.  The  recently  developed  parallel  propagation 
approximation  replaces  the  integration  over  wavenormal  distribution  with  a  closed  form 
expression,  and  can  be  quite  accurate.  However,  it  can  also  perform  badly,  especially 
for  electron  energy  >  1  MeV.  Here,  the  accuracy,  justification,  and  limits  of  the 
approximation  are  explored,  and  an  improved  version  is  presented.  It  is  based  on  a 
previously  developed  procedure  for  identifying  wavenormal  angles  compatible  with 
imposed  cutoffs  on  the  wave  frequency  distribution.  Because  it  also  requires  evaluation  at 
only  a  small  number  of  points,  it  features  computational  efficiency  comparable  to  the 
parallel  propagation  version,  while  preserving  contributions  from  oblique  waves  and  all 
relevant  harmonic  numbers.  Detailed  comparisons  are  presented  using  an  established 
model  of  nightside  chorus. 
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1.  Introduction 

[2]  Quasi-linear  diffusion  by  cyclotron-resonant  plasma 
waves  [Horne  and  Thorne ,  1998]  have  been  established 
as  a  likely  key  ingredient  of  the  behavior  of  electrons  in 
the  Earth’s  radiation  belts.  In  particular,  attention  has 
focused  on  whistler  mode  chorus  [Horne  et  al. ,  2005], 
hiss  [Meredith  et  al .,  2006],  and  electromagnetic  ion 
cyclotron  (EMIC)  waves  [Summers  and  Thorne ,  2003], 
alone  and  in  combination  [Summers  et  al .,  1998].  In 
conjunction  with  increasingly  compelling  empirical  stud¬ 
ies  [e.g.,  Meredith  et  al .,  2003a,  2003b;  lies  et  al .,  2006], 
multidimensional  diffusion  simulations  are  being  devel¬ 
oped  [Albert  and  Young ,  2005;  Varotsou  et  al .,  2005; 
Shprits  et  al.,  2006a].  Such  simulations  require  diffusion 
coefficients  evaluated  for  many  values  of  energy,  equato¬ 
rial  pitch  angle  a0,  and  (in  3D)  L  value,  under  a  wide 
variety  of  evolving  geophysical  conditions.  Furthermore, 
the  simulation  of,  say,  1  MeV  electrons  requires  a 
computational  domain  (and  diffusion  coefficients)  extend¬ 
ing  significantly  beyond  1  MeV,  where  boundary  con¬ 
ditions  (usually  /  =  0)  are  applied.  Thus  the  diffusion 
coefficients  must  be  evaluated  reasonably  quickly,  as  well 
as  reasonably  accurately. 

[3]  Toward  this  end.  Summers  [2005,  2007a,  2007b] 
have  developed  the  parallel  propagation  approximation, 
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which  collapses  the  usual  broad  wavenormal  distribution 
of  the  waves  to  a  magnetic  field-aligned  (or  anti-aligned) 
population.  The  integration  over  wavenormal  angle  0  of 
quasi-linear  theory,  which  picks  the  resonant  waves  out  of 
the  distribution,  is  replaced  by  a  closed  form  expression 
which  is,  naturally,  much  faster  and  easier  to  evaluate. 
The  neglect  of  oblique  waves  also  removes  contributions 
with  resonant  harmonic  number  n  other  than  -1.  Con¬ 
trary  to  long-established  expectations  based  on  plasma- 
spheric  hiss  [Lyons  et  al.,  1972;  Albert ,  1994,  1999], 
preliminary  comparison  of  this  approach  with  full  quasi- 
linear  calculations  [Shprits  et  al.,  2006b]  indicated  sur¬ 
prisingly  good  agreement,  with  one  caveat  discussed  in 
section  2.4. 

[4]  This  paper  explores  further  the  accuracy  of  the  parallel 
propagation  approximation,  interprets  it  as  a  weighted  aver¬ 
age  over  wavenormal  angle  9  of  the  full  expression,  shows 
where  it  is  likely  to  perform  poorly,  and  suggests  an  improved 
version.  This  modification  is  based  on  previously  developed 
procedures  for  identifying  wavenormal  angles  compatible 
with  frequency  cutoffs  imposed  on  the  modeled  wave  distri¬ 
bution.  It  also  replaces  an  integration  over  wavenormal  angle 
with  evaluation  at  a  small  number  of  values,  while  retaining 
contributions  from  oblique  waves  and  relevant  harmonic 
resonance  numbers.  The  contribution  from  n  =  —  1  is  as  good 
an  approximation  to  the  true  value  as  the  parallel  propagation 
result  is,  while  the  additional  terms  estimate  their  counter¬ 
parts  in  the  full  calculation.  The  result  is  a  more  favorable 
compromise  between  computational  speed  and  qualitative 
accuracy.  These  features  are  demonstrated  for  several  ver¬ 
sions  of  calculations  using  an  established  model  of  nightside 
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chorus  appropriate  for  L  —  4.5,  for  a  wide  range  of  energy  and 
equatorial  pitch  angle. 

2.  Diffusion  Coefficients 

[5]  With  slight  change  of  notation  from  Albert  [2005],  the 
local  pitch  angle  diffusion  coefficient  Dna  of  Lyons 
[1974b],  with  dimensions  of  p2/t,  can  be  written  as 

D<'“  _  ^wave  ry, 

p2  ^2  ft!  Z^  oa* 

*  1  o  11=—^ 


r  r/rnix 

££0  =  /  s\r\  6d6  A„  G\G2, 
Jomm 

with 


[8]  In  equation  ( 1 ),  G2(u;,  0)  and  A„(u ;,  0)  are  evaluated  at 
the  resonant  frequency  u  corresponding  to  0  and  ny  deter¬ 
mined  from  the  cyclotron  resonance  condition.  There  may 
be  several  such  values  of  a ;,  hence  the  sum  over  u  in  ( 1 ). 
Thus  for  each  ny  A„(u;,  0)  and  Gj(u;)  are  just  functions  of  0. 
However,  in  evaluating  A f(u>),  the  function  T  is  evaluated  at 
(u<0),  O').  The  functional  form  of  T  arises  from  the  Jacobian 
of  a  transformation  from  (kj_,  £||)  to  (a;,  0),  as  discussed  in 
the  Appendix  of  Lyons  [1974b]. 

[9]  The  purpose  of  writing  the  diffusion  coefficients  in 
the  form  equation  (1)  is  to  try  to  isolate  the  dependence  of 
Daa  on  the  parameters  of  the  frequency  and  wavenormal 
distributions.  Thus  A„  and  T  are  independent  of  the 
modeled  wavefield,  G|  depends  only  on  the  frequency 
parameters,  and  (for  lu  within  the  frequency  cutoffs)  G2 
and  N(uj)  depend  only  on  the  wavenormal  parameters. 


7r  sec  9  (~  sin2  a  +  snQc/urf)2 

2|v„/c|3  "  \\-{du/dk^)0/n\ 


(2) 


2.1.  Resonances 

[10]  The  resonance  condition  is  usually  given  as 


a;  -  *||V||  =  snQc/r 


(5) 


Gi  M 


qcb2{u) 


G2  (cj,  0)  = 


gM 


dO'  sin  0Tg>i..(0'), 


r  =  p 


Op 

duj 


However,  it  is  also  usually  assumed  that  the  wave  intensity 

(3)  is  distributed  symmetrically  with  respect  to  the  sign  of  k\\9  so 
that  the  particle  encounters  equatorward  and  anti-equator- 
ward  waves  at  eveiy  instant  [Lyons  et  al.y  1971].  This  is 
automatically  accounted  for  by  solving 

(4) 

(u-snQJ  7)2=  (Vn)~.  (6) 


For  Dnp  (or  Dpp)y  the  only  change  is  that  A„(u;,  0)  is 
multiplied  by  one  (or  two)  additional  factors  of  [sinacosa/ 
(— sin~a  +  snQcJury)]. 

[6]  In  these  formulas  p  is  the  refractive  index,  kc/cjy  given 
as  a  function  of  (uy  0)  by  cold  plasma  theory  [e.g.,  Stix , 
1992].  The  term  <P;,  accounts  for  the  relationships  between 
the  components  of  wave  electric  and  magnetic  fields  as  well 
as  the  gyroaveraging  of  the  resonant  wave-particle  phase, 
and  is  given  by  equation  (9)  of  Lyons  [1974b];  it  contains 
Bessel  functions  Jn  and  Jn±{  with  argument  5  k±  p±/mQc 
(where  s  is  the  sign  of  the  particle’s  charge,  m  is  its  rest 
mass,  7  is  its  relativistic  factor,  and  Qc  is  the  absolute  value 
of  its  nonrelativistic  cyclotron  frequency  in  the  background 
magnetic  field  B0).  Mathematically,  the  Bessel  functions  are 
appreciable  (and  oscillatory)  when  the  argument  is  >|w|. 
Since  the  argument  can  be  written  as  (sury/Q.c  -  w)tana 
tan0,  contributions  from  large  |/i|  are  obviously  favored  by 
large  0,  large  energy,  and  resonant  waves  at  the  latitude  of 
the  mirror  point  (where  a  is  maximum).  In  particular, 
and  A„  are  zero  at  0  =  0  unless  n  =  -1. 

[7]  The  function  B2(uj)  describes  the  frequency  distribu¬ 
tion  of  wave  power,  and  is  taken  to  be  zero  unless  uj  lies 
between  the  lower  and  upper  cutoffs  wLC  and  ljuc.  Simi¬ 
larly,  the  distribution  of  wave  power  with  wavenormal  angle 
0  is  described  by  £^(0),  which  is  zero  unless  0  lies  between 
^min  0max-  If  is  standard,  though  not  necessary,  to  use 
truncated  Gaussian  distributions:  B2(uj)  =  exp[-(u;  - 
u>,n)2/6uf]  and  gJ9)  =  exp[-(tan<?  -  tan0m)2/tan2(<56>)]. 
The  effect  of  using  a  Gaussian  in  tan#  instead  of  0  in  the 
model  is  to  reduce  the  power  ascribed  to  large  values  of  0. 


rather  than  equation  (5),  for  frequencies  of  resonant  waves 
in  either  direction.  These  resonances  all  occur  in  1/4 
bounce,  from  the  equator  to  the  mirror  point,  without 
accounting  for  the  sign  of  *„V||  [Lyons  et  al.,  1972].  Note 
that  the  integrals  in  equation  (1)  and  equation  (4)  are  only 
over  positive  values  of  0,  which  might  more  properly  be 
written  as  |0|.  If,  instead,  the  same  total  wave  intensity  is 
distributed  among  only  anti-equatorward  waves  (as  is 
probably  more  realistic  for  chorus  and  EMIC  waves),  the 
particle  takes  twice  as  long  (going  from  equator  to  mirror 
point  and  back)  to  encounter  resonances  with  either  sign  of 
£||V||.  However,  for  the  same  total  intensity  Bi avc»  each 
resonant  wave  will  have  twice  the  power  as  in  the 
bidirectional  case,  so  the  bounce  averaged  result  will  be 
the  same.  The  prescription  for  bounce  averaging,  while 
transforming  from  local  a  to  equatorial  ao,  was  given  by 
Lyons  et  al.  [1972]. 

2.2.  Behavior  of  the  Terms 

[n]  The  behavior  of  the  various  terms  vs.  0  is  shown  in 
Figures  1  and  2,  for  1  MeV  electrons  with  a  =  60°  and 
UpJ£le  =  2.5  (given  this  ratio,  it  is  not  necessary  to  specify  the 
L  value).  The  top  panels  show  A/f(o;,  0)  and  also  As 
illustrated  for  both  n  =  -1  and  n  =  -2,  the  Bessel  function 
term  $1  tends  to  be  the  controlling  factor  of  A„. 

[12]  The  middle  panels  show  u<0)  and  also  G,(cj).  To 
calculate  Gi,  the  frequency  parameters  were  taken  to  be 
~  0.35f2e,  Slj  =  0.15f2e,  ujlc  =  0.05Qe,  and  uj^jc  = 
0.65f2e,  as  in  the  nightside  chorus  model  of  Horne  et  al. 
[2003].  Note  that  the  variation  of  u  with  0  is  suprisingly 
weak  for  small  and  moderate  0,  and  lies  well  within  the 
frequency  cutoffs.  This  may  be  rationalized  by  considering 
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Figure  1.  The  behavior  of  various  factors  of  the  diffusion 
coefficients  for  1  MeV  electrons,  with  a  =  60°,  wpe/Qe  =  2.5 
and  n  =  —  1.  The  various  frequency  and  wavenormal 
parameters,  for  nightside  chorus,  are  given  in  the  text.  Top 
panel:  A „(uj,  0)  (left  vertical  axis)  and  <£,?  (right  vertical 
axis).  Middle  panel:  i j(0)  (left  axis)  and  G|(u;)  (right  axis). 
Bottom  panel:  G2  (left  axis)  and  G|G2  A„  sin#  (right  axis). 


the  resonance  condition  in  the  form  V(cj,  0)  =  ^(cj,  # ),  where 
^  =  u r/k2 c2,  as  discussed  in  detail  by  Albert  [2005].  Plotted 
against  lu  at  fixed  # ,  the  curves  of  Fand  'I'  both  decrease  as  # 
increases,  so  that  the  location  of  their  intersection  tends  not 
to  change  much.  From  another  point  of  view,  the  simplified 
treatment  of  Lyons  et  al.  [1972]  gave  u ;  ~  sec 20  (or  uj  ~ 
cos 2 6  for  n  =  0),  so  that  uj  is  nearly  constant  (du/dO  is  small) 
for  small  #.  This  near  constancy  of  lj(0)  will  be  exploited  in 
the  next  section. 

[13]  The  bottom  panels  of  Figures  1  and  2  show  G2  and 
the  product  G|G2A„  sin#,  whose  integral  is  D£tt.  The  wave- 
normal  parameters  #„,  =  0 ,60  =  30°,  #min  =  0, 0max  =  45°  were 
used,  also  as  done  by  Horne  et  al.  [2003].  It  is  seen  that  the 
interplay  between  4>„  and  G2  sin#  controls  which  6  values 
contribute  to 

2.3.  Interpretation  as  a  6  Average 

[14]  As  noted  above,  the  middle  panels  of  Figures  1  and  2 
show  that  for  small  and  moderate  6  values  lu(0)  is  roughly 
constant.  Thus  in  evaluating  the  normalization  N( u),  it  may 
be  acceptable  to  replace  r(u<0).  O')  with  TM#'),  O').  This 


approximation  has  the  benefit  that  N  can  then  be  taken 
outside  the  0  integral  in  equation  (1),  giving 


j£r  ddsm0&„GlgJ6) 

JtZ  sintfr fcW 


(7) 


This  form  may  be  interpreted  as  a  weighted  average: 


a: 


(8) 


with 


f'r  d6 sin  9r(uj(6) ,  0)g^(8)f(6) 

<f(0))  =  - -  —  -  .  (9) 

dosmor(u>(o).e)gj0) 

Figure  3  shows  A„  G\/T  and  the  weighting  function  sin# 
Fg^(#)  corresponding  to  Figures  1  and  2.  Note  that  for  n  = 
—  1,  A„  G,/r  is  practically  constant  for  small  #,  while 
much  of  the  #  dependence  in  equation  (1)  has  been 
isolated  in  the  weighting  function.  Thus  the  weighted 
average  of  the  relatively  constant  factors  is  expected  to  be 
insensitive  to  the  exact  form  of  the  weighting  ftinction, 
and  to  the  #  range  over  which  the  average  is  performed. 
The  actual  values  of  the  weighted  average  are  indicated  by 
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Figure  2.  Same  as  Figure  1,  for  n  =  -2. 
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Figure  3.  The  combination  A„(7|/r  and  weighting  func¬ 
tion  sin0  1^,(0)  corresponding  to  Figures  1  and  2.  The  local 
diffusion  coefficient  Dnan  is  approximated  by  the  integral 
over  0  of  their  product  in  equation  (8).  The  weighted 
average  of  A„G|/r  is  indicated  by  the  dotted  curves,  and  its 
value  evaluated  at  0,.  is  shown  by  the  triangles. 


the  dotted  lines,  while  the  points  marked  by  triangles  will 
be  discussed  below. 

2.4.  Evaluation  at  Discrete  0  Values 

[15]  The  well-known  mean  value  theorem  [e.g.,  Gradshteyn 
and  Ryzhik ,  1980]  guarantees  that  the  average  of  the  expres¬ 
sion  in  equation  (8)  can  be  obtained  by  simply  evaluating  it  at  a 
single  (unknown)  value  0O  between  0min  and  0max: 


00) 

With  0O  =  0,  this  result  and  the  corresponding  ones  for  D«p 
and  Dnpp  are  exactly  1/2  times  those  of  equations  (5-7)  of 
Summers  et  al.  [2007a]. 

[i6]  Full  numerical  integration  with  the  PADIE  code 
[Glauert  and  Horne ,  2005]  also  manifested  this  numerical 
factor  relative  to  Summers  [R.  Home,  personal  communi¬ 
cation,  2006],  which  is  evident,  though  not  noted,  in  the 
plots  of  Shprits  et  al.  [2006b].  The  descrepancy  is  tenta¬ 
tively  ascribed  to  the  explicit  factor  of  2  in  equation  1  of 
Summers  [2005],  which  is  not  apparent  in  earlier  work 


[Kennel  and  Engelmann ,  1966,  equation  4.1;  Lerche , 
1968,  equation  14;  Lyons  et  al .,  1971,  equation  1;  Lyons, 
1974a,  equation  1],  and  not  to  any  error  in  the  ensuing 
derivation.  Pending  the  ultimate  resolution  of  this  issue,  this 
paper  will  consider  the  parallel  propagation  approximation 
in  the  form  of  equation  (10),  with  0O  =  0. 

[i  ?]  More  generally,  if  equation  (10)  is  to  be  a  reasonable 
approximation  to  equation  (8),  the  condition  ujlc  <  u<0o)  < 
iOyC  is  required;  values  of  uj(0)  which  satisfy  the  resonance 
condition  but  do  not  lie  within  the  cutoffs  are  termed 
noncontributing.  Techniques  for  identifying  such  values  of 
0  have  previously  been  developed  for  whistler  mode  waves 
[Albert,  1999,  2004,  2005],  EMIC  waves  [Albert,  2003], 
and  superluminous  and  Z  mode  waves  [Albert,  2007],  for 
the  purpose  of  performing  the  full  integral  equation  (1) 
efficiently,  and  are  reviewed  in  Appendix  A.  These  techni¬ 
ques  yield  0  intervals,  or  ranges,  (0a,  0h),  which  can  be 
directly  applied  to  choosing  where  to  estimate  equation  (8). 

[is]  If  there  is  only  one  0  range,  0O  is  taken  as  some 
representative  value  such  as  (0a  +  0h)/2  or  tan_,[(tan0rl  + 
tan0,,)/2]  (as  used  below  for  convenience).  If  there  is  more 
than  one  0  range,  (7-9)  are  expressed  as  sums  of  averages 
over  each  range  r: 


ry,  ^  d0  sin  0A„G|g^.(0) 

frd0  sin  erg  JO) 


/  A„G|  \  ^  A„G| 


(id 


where 


a  _  f^0sm0TgJfi)  _  jrJ6sm6gJ0) 

'  Er  frd0sinOrgu(O)  ~  Zr  l desmOg J0)  (l2) 

This  few-point  estimate  captures  at  least  some  measure  of 
the  contribution  to  equation  (1)  by  oblique  waves  and  all  n 
values.  In  favorable  situations  the  0  ranges  will  be  narrow, 
so  that  A„G,/T  is  nearly  constant  within  each  one;  at  worst, 
the  procedure  defaults  to  0r  =  (0min  +  0max)/2.  Since  each 
range  may  contain  some  noncontributing  values  as  well  as 
all  contributing  values,  it  is  possible  for  the  chosen  0,.  to 
misrepresent  an  entire  range  as  noncontributing;  however, 
the  techniques  discussed  in  Appendix  A  are  discriminating 
enough  that  this  is  rare.  In  Figure  3,  the  values  of  A„G,/r 
evaluated  at  0r  are  marked  by  triangles  and  give  good 
approximations  to  the  weighted  averages,  indicated  by  the 
dotted  lines. 

[19]  The  fractions  ar  in  equation  (12)  may  be  evaluated 
numerically.  However,  with  gjfi)  given  by  the  truncated 
gaussian  of  x  =  tan0  mentioned  above,  it  is  convenient  to 
further  approximate  sin0  dO  *  xdx/(  1  +  x2)3'2  by  simply  xdx. 
For  x  <  1  this  does  little  harm,  as  shown  in  Figure  4,  and 
allows  ar  to  be  evaluated  analytically.  This  procedure  was 
followed  in  the  numerical  evaluations  below. 


3.  Results 

[20]  Figure  5  shows  several  versions  of  bounce-averaged 
diffusion  coefficients  calculated  with  u )pj£le  =  2.5,  ujm  = 

0.350,  6u  =  0.150,  vlc  =  0.050,  Vue  =  0.650,  0m  =  0, 
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n  =  -  1 


to  the  full  results  (black)  for  0.3  MeV,  but  are  only  fair 
approximations  to  the  full  result  for  n  =  -1  (green)  for 
higher  energy  and  a0  above  60°,  and  poor  approximations 
to  the  full  results.  Finally,  the  magenta  curve  shows  equa¬ 
tion  (11)  again  evaluated  at  0r  but  summed  over  n  =  -5  to  5, 
which  seems  to  provide  a  good  approximation  to  the  full 
results  (black). 

[22]  Figure  6  shows  the  same  four  approximations,  nor¬ 
malized  by  the  full  calculation.  Again,  the  green  curves 
show  equation  (1)  with  n  =  —  1  only,  the  red  curves  show 
equation  (10)  with  #0  =  0,  and  the  blue  curves  show 
equation  (11)  with  values  of  0r  but  n  =  -1  only.  The 
magenta  curve,  showing  equation  (11)  with  the  calculated 
values  of  0r  and  w  =  —  5  to  5,  clearly  gives  the  best 
approximation  (to  a  ratio  of  1 ,  shown  by  the  dashed  black 
line),  with  the  exception  of  the  region  a0  >  75°,  E  =  1  MeV. 
Note  that  the  region  of  good  agreement  in  Dpp  for  the  red 
curve  at  1  MeV  and  3  MeV  between  cv0  =  60°  and  75°  is 


e 

Figure  4.  The  weighting  function  sin#  Tg^f#)  of  Figure  3 
(solid  curve),  and  the  approximations  sin#g^(#)  (dashed 
curve)  and  sintfgj^/cos3#  (dotted  curve),  where  gu(0)  is  a 
Gaussian  function  of  tan#.  The  two  approximations  have 
been  scaled  to  have  the  same  maximum  as  the  exact 
expression;  the  scale  factor  cancels  in  equation  (12).  The 
second  approximation  allows  (12)  to  be  evaluated 
analytically. 


60  =  30°,  #min  =  0,  #max  =  45°  as  above,  £wavc  =  50  pT,  and 
latitude  range  0  <  A  <  15°  [Horne  et  ai ,  2005].  The  black 
curves  show  the  full  calculation  of  equation  (1),  including 
all  n  between  —5  and  5,  which  will  be  regarded  as  “the  right 
answer,”  to  be  approximated. 

[21]  The  green  curves  show  equation  (1)  calculated  with 
just  n  =  -1,  and  are  very  close  to  the  full  results  for 
0.3  MeV,  but  miss  dominant  contributions  from  other  n 
values  for  1  MeV  and  especially  3  MeV  below  a0  =  60°. 
The  red  curves  show  equation  (10)  evaluated  at  #0  =  0 
(hence  n  =  —  1  only),  as  in  Shprits  et  al.  [2006b],  while  the 
blue  curves  show  equation  (11)  evaluated  with  tan#r  at  the 
center  of  the  computed  #  ranges,  but  only  for  n  =  —  1 .  They 
(red  and  blue)  are  always  close  to  each  other  and  very  close 
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Figure  5.  Several  versions  of  bounce-averaged  diffusion 
coefficients  (in  sec-1)  for  0.3  MeV,  1  MeV,  and  3  MeV 
electrons  and  nightside  chorus,  with  wave  parameters  given 
in  the  text.  The  black  curves  show  the  full  calculation  with 
n  <5.  The  green  curves  show  the  full  calculation  with  just 
n  =  —  1.  The  red  curves  show  results  using  the  parallel 
propagation  approximation.  The  blue  curves  show  the  “# 
range”  approximation  discussed  in  the  text,  but  only  with 
n  =  -1.  The  magenta  curves  show  the  “#  range” 
approximation  with  n  <  5. 
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Figure  6.  The  four  approximations  shown  in  Figure  5, 
normalized  by  the  results  of  the  full  calculation. 


fortuitous  because  the  green  curve,  computed  from 
equation  (1)  with  n  =  —  1  only,  also  differs  greatly,  so  that 
the  diffusion  actually  comes  from  other  n  values,  which  the 


red  (and  blue)  curves  do  not  include.  Also,  the  corresponding 
agreement  in  Dftoao  is  poor  for  all  but  the  magenta  curve. 

3.1.  Sensitivity  to  Parameter  Changes 

[23]  Following  Shprits  et  al.  [2006b],  Figure  7  shows  the 
result  of  doubling  the  density  ratio  from  uope/Qe  =  2.5  (solid 
curves)  to  5.0  (dashed  curves).  The  full  calculations  for 
1  MeV,  shown  in  black,  show  that  the  diffusion  coefficients 
decrease  at  every  a0.  The  calculations  using  6r  (blue) 
decrease  as  well,  maintaining  the  same,  mostly  good,  degree 
of  approximation.  The  calculations  using  0o  =  0  (red),  which 
were  only  fair  to  begin  with,  become  worse  at  the  higher 
density. 

[24]  Figure  8  shows  the  effect  of  increasing  the  latitudinal 
extent  of  the  waves,  from  15°  (solid  curves)  to  25°  (dashed 
curves).  None  of  the  calculations  are  affected  for  a0  >  60°, 
since  these  particles  mirror  below  15°  latitude.  For  smaller 
a0,  the  calculations  with  0O  =  0  (red)  improve  their 
agreement  with  the  full  calculations  (black),  but  the  calcu¬ 
lations  using  0r  (blue)  improve  even  more,  becoming  almost 
indistinguishable  from  the  full  results  with  the  wider  latitude 
range. 

[25]  Shprits  et  al.  [2006b]  computed  diffusion  coeffi¬ 
cients  for  different  parameter  values  with  90  =  0  and 
compared  them  to  each  other,  and  suggested  that  the 
variation  of  the  results  with  the  (imperfectly  known)  pa¬ 
rameter  values  outweighed  the  benefit  of  performing  the  full 
calculations.  However,  this  is  not  conclusive,  because 
comparison  was  not  made  to  full  calculations  using  the 
same  parameter  values.  The  present  results  show  that  the 
parallel  propagation  approximation  results  are  significantly 
more  sensitive  to  the  parameter  changes  than  the  full  results 
are,  while  the  approximations  based  on  the  0  ranges 
reproduce  the  behavior  of  the  full  calculations  with  good 
fidelity. 

3.2.  Qualitative  Comparison  at  L  =  4.5 

[26]  Figures  9-11  show  diffusion  coefficients  calculated 
with  the  same  wave  model  but  with  tJpe/Qe  =  3.43  on  the 


Figure  7.  The  effect  of  doubling  the  density  ratio,  from 
upJSle  =  2.5  (solid  curves)  to  5.0  (dashed  curves),  on 
diffusion  rates  (in  sec-1)  for  1  MeV  electrons.  The  full 
calculations  are  shown  in  black,  the  parallel  propagation 
approximation  in  red,  and  the  “0  range”  approximation  in 
blue. 


Figure  8.  The  effect  of  increasing  the  maximum  latitu¬ 
dinal  extent  of  nightside  chorus  from  15°  (solid  curves)  to 
25°  (dashed  curves).  The  full  calculations  are  shown  in 
black,  the  parallel  propagation  approximation  in  red,  and  the 
“0  range”  approximation  in  blue.  (The  blue  and  black 
dashed  curves  are  virtually  indistinguishable.) 
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Figure  9.  Diffusion  coefficients  (in  s~‘)  for  electrons  and 
nightside  chorus  waves  with  cupe/Qe  =  3.43,  for  £=0.1  MeV 
to  5  MeV,  from  a  full  calculation  with  n  <  5. 

equator,  as  appropriate  for  L  =  4.5  [Horne  et  al .,  2005].  The 
calculations  were  done  for  85  values  of  a0  from  5°  (just 
outside  the  loss  cone)  to  89°,  and  49  logarithmically  spaced 
values  of  energy  between  0. 1  MeV  and  5  MeV  (inclusive). 
Figure  9  shows  the  full  calculations  with  n  between  —5  and 
5,  which  took  about  5  hours  of  computer  time.  Dark  blue 


regions  at  the  upper  left  comer  (large  £,  small  ao)  of  panels 
(a),  (b),  and  (c)  are  off  the  color  scale,  but  are  not  zero. 
Figure  10  shows  the  results  using  the  approximation  equa¬ 
tion  (10)  with  at  0O  =  0  (and,  necessarily,  n  =  —  1  only). 
They  took  only  about  4  seconds  to  obtain,  and  show  broad 
areas  of  qualitative  agreement  with  Figure  9,  but  also  show 
large  areas  of  disagreement.  The  uncolored  areas  in  the 
upper  left  of  all  four  panels  indicate  diffusion  coefficients  of 
identically  zero,  indicating  that  no  (n  =  -1)  resonances 
were  present.  These  regions  include,  but  extend  beyond,  the 
dark  blue  small-value  regions  of  Figure  9.  Also,  the  cross 
diffusion  coefficient  is  often  estimated  poorly,  in  that  the 
boundary  separating  positive  from  negative  values  does  not 
match  Figure  10  at  all  well.  Finally,  Figure  11  shows  the 
results  using  equation  (11)  using  6n  for  n  =  -5  to  5.  These 
calculations  took  about  90  seconds,  and  show  much  better 
agreement  with  Figure  9  than  Figure  10  does.  Finite  values 
are  reported  at  all  values  of(£,  a0),  and  the  general  behavior 
of  DttoP  is  captured  well. 

4.  Conclusions 

[27]  The  local  quasi-linear  diffusion  coefficients  have  been 
approximated  as  weighted  averages  over  wavenormal  angle 
6.  This  relies  on  the  resonant  frequency  uj  being  nearly 
constant  as  a  function  of  wavenormal  angle  0.  Naturally,  this 
is  favored  when  0max  is  close  to  0min  (so  0max  should  be  small 
if  0min  is  small). 

[28]  Using  previously  developed  techniques,  the  sub¬ 
ranges  of  0  which  contribute  to  the  diffusion  can  be 
identified  (see  Appendix  A),  and  the  averages  over  each 
range  can  be  further  approximated  by  evaluation  of 
equation  (11)  at  single  values  9n  This  retains  the  effect  of 
oblique  waves,  and  their  resonances  with  arbitrary  harmonic 
number  n.  The  parallel  propagation  formulas  [ Summers , 
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Figure  10.  Same  as  Figure  9,  but  using  the  parallel 
propagation  approximation. 


Figure  11.  Same  as  Figure  9,  but  using  the  “0  range” 
approximation  with  n  <  5. 
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2005,  2007a]  are  recovered  by  evaluation  at  6  =  0  (hence 
n  =  ~  1  only),  except  for  a  factor  of  2  as  discussed.  This 
approach  may  also  be  applied  to  EMIC  waves  [Albert, 
2003]  as  well  as  to  superluminous  (L-X,  L-O,  and  R-X) 
and  Z  mode  waves  which  can  be  quite  oblique  [Summers  et 
al.,  2001;  Albert ,  2007],  though  this  is  deferred  to  future 
work. 

[29]  The  full  integration  of  equation  (1),  as  in  Albert 
[2005]  or  Horne  et  al.  [2005],  remains  the  “gold  standard” 
for  such  calculations,  though  they  are  moderately  demand¬ 
ing  of  computer  time.  The  parallel  propagation  expressions 
are  much  faster  and  easier  to  evaluate,  and  for  typical  chorus 
wave  parameters  give  reasonable  estimates  of  the  n  =  -1 
contribution,  when  present.  However,  this  is  not  always  the 
dominant  term,  especially  for  E  >  1  MeV.  The  scheme 
presented  here  may  be  considered  a  hybrid  of  the 
approaches  of  Summers  [2005]  and  Albert  [2005].  It  retains 
much  of  the  convenience  and  most  of  the  speed  of  the 
parallel  propagation  approximation  with  greater  accuracy,  as 
demonstrated  for  >  1  MeV  electrons  interacting  with  chorus. 
Thus  it  may  be  an  optimal  approach  for  practical  but 
comprehensive  simulation  of  quasilinear  diffusion  in  the 
radiation  belts. 

Appendix  A:  Finding  Resonant  0  Ranges 

[30]  To  find  conditions  for  resonant  6  ranges  consistent 
with  lower  and  upper  frequency  cutoffs,  the  resonance 
condition  equation  (6)  is  written  in  some  form  6)  = 
g(^>  0),  and  graphs  of / and  g  are  analyzed  geometrically  as 
functions  of  u  at  fixed  0.  Within  the  interval  ( uLC ,  u joc),  the 
curves  o f/(u/)  and  g(u/)  cannot  intersect  unless  (i)  max(/)  > 
min(g)  (otherwise /always  lies  below  g ),  and  (ii)  min(/)  < 
max(g)  (otherwise  f  always  lies  above  g).  Imposing  con¬ 
ditions  (i)  and  (ii)  retains  all  0  values  for  which  u/(0)  lies 
within  the  frequency  cutoffs,  while  eliminating  many 
though  not  necessarily  all  other  0  values.  Of  course,  the 
6  values  are  also  required  to  lie  between  0min  and  0max. 

[31]  The  analysis  has  previously  been  performed  for  the 
resonance  condition  in  the  form 


k2c2  c2  (  \2  cos2  UJn  ~  snQch  (Al) 

K  C  (u-  LUn) 

[Albert,  1999,  2003,  2005,  2007],  because  of  the  relatively 
simple  geometric  behavior  of  these  expressions  as  functions 
of  u.  The  resonance  condition  can  also  be  written  as 

*V  0)  =  (u-  u,,)2  ^  =  IV(U,  6).  (A2) 

1 

This  form  is  even  more  convenient  because  k1  is  a 
monotonically  increasing  function  of  u  at  fixed  6,  for  all 
cold  plasma  wave  modes,  while  fV(uj)  also  has  simple 
behavior  [Albert,  2007].  Thus  for  a <  0,  condition  (i) 
above  becomes  >  W(lolc)>  while  condition  (ii) 

becomes  <  (us uc )•  For  us„  >0,  there  are  three 

subcases.  When  u /„  <  ujlc,  the  conditions  are  the  same  as  just 
given.  When  u j„  >  a >uc,  the  conditions  become  ^(usuc)  > 
W(louc)  and  k2(usLC)  <  W(ijjLC ).  Finally,  when  ujlc  <  us„  <  usuc , 
usuc*  condition  (i)  is  automatically  satisfied  since  min(Jf')  = 


0,  while  condition  (ii)  requires  either  k2(uiLC)  <  W(usLC)  or 
k~(^Lc)  <  **W).  Each  of  these  inequalities  is  easily 
converted  to  quadratic  form,  schematically 

A  cos4  0  4-  B cos2  0  +  C  >  0  (A3) 

(where  the  coefficients  depend  on  usLC  and  ljuc ),  so  that  0 
ranges  of  the  form  0a  <  0  <  0h  are  easily  calculated  from  the 
roots  of  the  quadratic.  Limited  numerical  experiments 
indicate  that  these  conditions  based  on  (A2)  are  about  as 
effective  as  corresponding  conditions  based  on  (Al),  while 
being  simpler  to  implement.  Since  they  are  not  equivalent, 
both  may  be  used. 
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